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ABSTRACT 

This paper describes a numerical scheme for muhi-fluid hydrodynamics in the hmit 
of small mass densities of the charged particles. The inertia of the charged particles 
can then be neglected, which makes it possible to write an evolution equation for 
the magnetic field that can be solved using an implicit scheme. This avoids the severe 
restriction on the stable timestep that would otherwise arise at high resolution, or when 
the Hall effect is large. Numerical tests show that the scheme can accurately model 
steady multi-fluid shock structures both with and without sub-shocks. Although the 
emphasis is on shocks in molecular clouds, a multi-dimensional version of this code 
could be applied to any Astrophysical flow in which ambi-polar diffusion or the Hall 
effect, or both play a significant role. 
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1 INTRODUCTION 

In dense molecular clouds, the density of charged particles 
can be so low that the scale on which ambi-polar resistivity 
becomes important cannot be assumed to be neglibly small. 
The most obvious effect of this is that it allows the existence 
of shock structures that are much thicker than those deter- 
mined by viscous effects (see e.g. Draine & McKee 1993), 
but the enhanced magnetic diffusion also plays a role in the 
large scale dynamics of such clouds (see e.g. Mouschovias 
1991). 

The structure of shocks in which the dissipation is 
due to ambi-polar resistivity rather than viscosity has been 
studied extensively (e.g. MuUan 1971; Draine 1980; Flower, 
Pineau des Forets & Hartquist 1985; ; Draine 1986; War- 
die & Draine 1987; Chernoff 1987; Roberge & Draine 1990; 
Pihpp & Hartquist 1994; Wardle 1998). There are several 
reasons for this, of which perhaps the most important is 
that the flow time through these shocks is long enough for 
radiation from molecules to maintain the gas at a low tem- 
perature even for shock speeds as high as 50 km s~^. In- 
deed, Wardle (1998) points out that this process is so effec- 
tive that such shocks are responsible for much of the infra- 
red H2 and CO line emission in molecular clouds (Draine & 
Roberge 1982; Chernoff, McKee & HoUerbach 1982; Smith 
& Brand 1990; Smith, Brand & Moorhouse 1991; Chrysos- 
tomou et al. 1997). The heating due to the currents also 
raises the temperature of the molecular gas to the point 
where a number of important chemical reactions proceed 
at a significant rate (Draine, Roberge & Dalgarno 1983; 
Flower, Pineau des Forets & Hartquist 1985; Pineau des 



Forets, Flower Hartquist & Dalgarno 1986; Draine & Katz 
1986a,b; Kaufman & Neufeld 1996a,b). 

Although the obvious way to determine the steady 
shock structures is to solve the steady equations directly, 
this is not a simple matter in the general case. There is no 
great difficulty when the steady solution can be obtained 
by integrating through the structure in the appropriate di- 
rection (e.g. Wardle & Draine 1987) and one can also deal 
with solutions containing sub-shocks if the system reduces 
to two differential equations (Chernoff 1987) . However, these 
are all special cases and the general case of oblique shocks 
with cooling and chemical reactions is a two point bound- 
ary value problem that can only be solved by a relaxation 
method such as that used by Draine (1980). In that case 
it is much simpler to use a time dependent code to obtain 
the steady solution. This also has the advantage that it can 
be used for unsteady problems in which these effects are 
important. 

A number of numerical schemes for the time dependent 
equations have been devised (Toth 1994; Smith & Mac Low 
1997; Mac Low & Smith 1997; Stone 1997; Chieze, Pineau 
de Forets & Flower 1998; Ciolek & Roberge 2002), but, with 
the exception of Ciolek & Roberge (2002), these all assume 
that all the charged particles have a large Hall parameter 
and are therefore tied to the magnetic field. The system 
can then be modeled as three fiuids: a neutral fluid, an ion 
fluid and an electron fluid. The ions and electrons may have 
different temperatures, but the magnetic field forces them 
to have the same velocity, which can, however, be different 
from that of the neutral fluid. 

In molecular clouds, the assumption that the Hall pa- 



© 0000 RAS 



2 S. A. E. G. Falle 



rameter is large is valid for electrons, ions and small grains 
such as polycyclic aromatic hydrocarbons (PAH), but not 
for grains whose radius, ttg, is larger than about 10~^ cm. 
Pilipp & Hartquist (1994) and Wardle (1998) have shown 
that a significant charge density of particles for which the 
Hall parameter is 0(1) has a profound effect on the shock 
structures. In particular, such particles induce a Hall resis- 
tivity which can lead to substantial rotation of the transverse 
field within the shock structure. For this reason, Ciolek & 
Robcrgc (2002) devised an algorithm for time dependent 
equations that include a significant charge density of parti- 
cles whose Hall parameter is not large enough for them to be 
tied to the field. However, as we shall see, the stable timcstep 
for their algorithm can become very small, especially if there 
is a significant Hall resistivity. This paper describes an al- 
gorithm that does not suffer from this restriction. 

Although the emphasis here is on shocks in molecular 
clouds, the algorithm is quite general and is efficient for 
any problem in which either ambi-polar diffusion or the Hall 
effect, or both are important. In particular, although ambi- 
polar diffusion in molecular clouds is due to drift between 
neutrals and charged particles, it can also arise due to drift 
between electrons and ions in plasmas with a small fraction 
of neutrals if the density is high enough for the ion Hall 
parameter to be 0(1) (e.g. Sano & Stone 2002). 



2 MULTI-FLUID EQUATIONS 

Consider a set of N fluids for which the generic one dimen- 
sional equations are [i = 1- ■ ■ N) 



dt dx 



^ + ^(p,«.q,+pa) 



dci d 



(1) 



aiPi(E + qi AB) 

(2) 



dt dx 



(3) 



+aipicii • (E -h Qi A B) 



Here — (ui,Vi,Wi) arc the velocities, a; the charge to 
mass ratios, ej the total energies and i the unit vector in the 
X direction. We have 

ei = Pi{Ui + iq?) 

where Ui is the internal energy per unit mass of fluid i. We 
shall assume that 



Ui = Wi + 



Pi 



where Wi is the energy associated with the internal states 
of the particles of fluid i. It can include ionization energy, 
chemical binding energy etc. 



The source terms are: 



mass transfer rate from fluid 



j to fluid i; Uj - momentum transfer rate from fluid j to fluid 
i; Gij - energy transfer rate from fluid j to fluid i; Hi - 
external energy source/sink for fluid i. Global mass, energy 



and momentum conservation clearly require that Sij = —sji, 
fij = —fji and Gij — —Gji. Draine (1986) derives expres- 
sions for these terms, but for our purpose it is sufflcient to 
consider their general form. 

The momentum transfer rate, fij, is 



fij — Cjj -|- SijClj 



Sji^ii 



where Cij describes collisions between the particles of flu- 
ids i and j. Since it involves binary collisions and tends to 
equalise the velocities of the two fluids, it must be of the 
form 

C,j = piPjK^j{Ti,Tj, |qj - qi|)(qj - q^), 

where Ti and Tj are the temperatures of fluids i and j. Since 
global momentum conservation requires Cij = —Cji, we 

must have K^, = K,i. 



The energy transfer rate, Gij, is 

— qi ' Cij -\- SijCj SjiCi ~\- Dij , 



where Dij describes the effect of collisions between the par- 
ticles of fluids i and j. Like Cij, it involves binary collisions 
and tends to produce equilibrium between the two fluids. It 
must therefore be of the form 

Dij = piPjLij{Ti,Tj, \qj - qi|). 

with Lij = when T, = Tj and q, = q^. Global energy 

conservation requires 



qi ■ Ci 



-G, 



The fields are determined from Maxwell's equations, 
which reduce to 



g-Bz _ dEy 

dt ' dt dx ' dt dx 
dBy _ dB^ _ 



dx 



dx 



(4) 



(5) 



since the displacement current may be neglected. The cur- 
rent, J, is given by 



J = ^ CtiPiCli. 



and we also have charge neutrality 



E 



CtiPi 



(6) 



(7) 



The units for E, B and J are such that the speed of light 

and the factor An do not appear. 

The above equations are extremely general and include 
those that other authors have used to model shocks in molec- 
ular clouds as special cases Despite this generality, they do 
assume that each fluid can be described by the hydrody- 
namic approximation, which is only true if the interactions 
between particles of the different fiuids are much weaker 
than those between particles of the same fluid. 

In the usual applications, the mass is dominated by a 
fluid, i = 1 say, consisting of neutral particles, which obvi- 
ously has zero charge to mass ratio. The other fluids are: an 
electron fluid, i = 2; an ion fluid, i = 3; fluids consisting of 
grains of various types and sizes, i = 4 - ■ ■ N. The validity 
of the hydrodynamic approximation for systems consisting 



© 0000 RAS, MNRAS 000, 000-000 



Multi-fluid Magnetohydrodynamics 3 



only of neutrals, electrons and ions has been discussed by 
Draino (1986) and Chernoff (1987). They conclude that the 
hydrodynamic approximation is valid for the neutrals and 
electrons, both of which should have a Maxwellian distribu- 
tion, but this may not be true for the ions. 

The grain fluids are obviously very different since they 
consist of particles with much larger masses than the other 
fluids. Their thermal velocity dispersion is therefore small 
compared to the drift velocities, which means that all the 
particles of a particular size and type have the same veloc- 
ity. The hydrodynamic approximation with zero pressure is 
therefore appropriate for the grain fluids. 

Although it is possible to construct a numerical scheme 
for these equations as they stand, this would not be suitable 
for the conditions in shocks in molecular clouds. Consider 
the neutral momentum equation and suppose that the drift 
velocities are small compared to the thermal velocity disper- 
sion of the neutrals, as they must be near the upstream or 
downstream ends of the shock structure. The collision term 
is then approximately 

JV JV 

^Cii = ^piPfA'H(Ti,Ti,0)(qi-qi). 

i=2 i=2 

If Qc is a typical velocity, such as the shock velocity 
relative to the upstream fluid, then a balance between this 
term and the inertia terms induces a length scale 

I ^ ^ <t 

It is clear that Id determines the shock thickness. Then since 
Pi ^ pi for all i > 1, the equivalent length scales for the 

other fluids arc 

Ic 



Ku{Ti,Tj,0)pr' 

which means that Id -C Id- 

The charged fluids are also acted on by electromagnetic 
forces, which induce a length scale 

I 

"~ aiB' 

which is the Larmor radius of particles with velocity Qc- 

Provided either Id ^ Id or lei ^ Ici, the inertia and 
pressure terms can be ignored in the momentum equation for 
fluid i. Since this is generally true, the momentum equations 
for all fluids except the neutral fluid reduce to 



aiPi{E + Hi A B) + f,i = 0, 



(8) 



where the fact that pi <^ pi for all i has allowed us to ignore 
momentum transfer from all but the neutral fluid. The same 
arguments tell us that in this case the energy equations for 
the charged fluids reduce to 



Hi + Gil + UiPiCii • (E qi A B) = 



(9) 



If lei Id, then particles i are closely tied to the field 
lines. It is convenient to define a Hall parameter, f5i, for each 
charged fluid by 



le 



UiB 



lei KliPl 

The ions and electrons have /3i ^ 1 and are therefore tied 
to the field, whereas the grains are not since they typically 
have /3i ~ 1. 



3 NUMERICAL METHOD 

Our method deals with the same equations as Ciolek & 
Roberge (2002), but uses a somewhat different approach, 
which, as we shall see, leads to a more robust and efficient 
scheme. 



3.1 Equations for the Neutral Fluid 

Since fij = —fji, the reduced momentum equations for the 
charged fluids, (8), the charge neutrality condition, (7) and 
the expression for the current, (6) give 



JV 



JV JV 



^ fli = - ^ fil = ^ a,p,(E q; A B) = J A B. 

i=2 i=2 i=2 

The momentum equation for the neutral fluid, (2), then be- 
comes 



dpiqi d 



+ (piwiqi + pii) = ^ fii = J A B. 



dt dx 



Similarly, the reduced energy equations, (9), give 



(10) 



JV 



JV 



^ Gli = - ^ Gil = ^ [Hi + aiPiCii ■ (E qi A B)] 



i=2 i=2 
JV 



= J-B + ^lHi + UiPiCii • (qi A B)]. 



<=2 



The energy equation for the neutrals, (3), then becomes 



dei d 



+ i: [^^i (ei +P^ + \p^<h)] = Hi+Y, Gii 



dt dx 



(11) 



: J • E ^ -h aip<q< • (qi A B)]. 



Equations (10), (11) and the continuity equation for the 
neutrals, (1) with i = 1, are just the ordinary gas dynamic 
equations with source terms. They can be written in the 
form 



^ + — =S 
dt dx 



(12) 



where 

Q = (pi,piqi,ei), 

is a vector of conserved variables, F is the vector of fluxes 
and S is the vector of source terms. We have F = F(Q), but 
the source term depends upon the magnetic fleld and the 
state of the other fluids. We therefore write 

S = S(Q,V) 

where V is a vector representing B, pi, q^ etc. 

3.2 Numerical Scheme for the Neutral Fluid 

Equations (12) are solved using the second order Godunov 
scheme described in Fallo (1991). This is a conservative 
scheme in which the numerical solution at time tn in the 
ceU with {j - 1/2) Ax < x < {j + 1/2) Ax is deflned to be 
the volume average 
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n 0+1/2) Ax 



\5{t„,x) dx. 



0-l/2)Ax 



Integrating (12) over the jth. cell and from t„ to 
t„ + At„, gives 

^(Qi - Q. ) + a;;:(f.+i/2 -^^-1/2) 



(13) 



1+1/2 



where ^^^1/2 are the time averages of the fluxes at the cell 

edges and S""''^''^ is the time average of the integral of the 
source term over the cell. 

We start by using equation (13) to compute a first order 
approximation, Q"^^^^, to the solution at the half time, 
^n+i/2 = in + Atn/2. In tMs step the fluxes and source term 
are approximated by 



5,n+l/2 

■ i+1/2 



F.(Q",Q"+l), S 



ri+1/2 



SCQ^V]"), 



where F,(QL,Qfl) is the flux in the resolved state for a 
gas dynamic Riemann problem for which the left and right 
states are Ql and Qr. 

To make the scheme second order in space, we use the 

solution at the half time to construct an average gradient of 
the primitive variables, P — (pi,qi,pi), in the ith cell 

/aP\"+^/^_ 1 , „+l/2 „+i/2 pn+1/2 pn+1/2-^ 



where the averaging function is 

r if a6 < 0, 

av{a,b) = < a^b + ab^ 

— o , .9 otherwise. 

The averaging function acts as a non-linoar switch that re- 
duces the scheme to first order in regions with large second 
derivatives, such as shocks etc (see e.g. van Leer 1977). 

The solution at = tn + At„ is then calculated from 
(13) with 

Fj+i/a =F*[Q(Pi),Q(P«)], 
where 



p ^ pn+l/2 



1, /9P\"+i/2 



and the source term given by 

gn+l/2 _ S(-Q"+l/2 Y»+l/2-j 

This is obviously not a complete scheme since we have 
yet to devise a method of advancing the variables, V, de- 
scribing the field and charged fluids. 



(Cowling 1957; Nakano & Umcbayaslii 1986). Here ro is the 
resistivity along the field, ri is the Hall resistivity and r2 
is the ambi-polar resistivity. If we define the conductivity 
parallel to the field, 

1 ^ 

i=2 

the Hall conductivity, 

N 

_ 1 Ciipi 

and the Pedersen conductivity, 

N 

BZ^(l + /3f) 
then the resistivities are 

ro = l/<To, n = fTi/(iTi + o-|), r2 = a-2/(<T? + cri) 

Substituting (14) into the Faraday equations (4) and 
using Ampere's law (5) to eliminate J gives the induction 
equation 



m dm 

dt dx 



— R— 

dx dx ' 



(15) 



where the hyperbolic flux is 
M = (0, uBy - vBx, uBz - wB^ 
and the resistance matrix is 
R = roRo + riRi + r2R2, 
with 



Ro 



Ri = 



R2 = 



/ 

V 

/ 

V 



S2 



B 



2 

V_ 

52 



B 



/ 1-^ 

B2 



v 



ByBz 

~B2~ 



ByBz 

£2 



£2 / 



Equation (15) can used to advance the field, but, as we 
shall see, an explicit approximation to this is very inefficient 
if the Hall term is much larger than the ambi-polar diffusion 
term. 



3.3 Equations for the Charged Fluids and Fields 

The reduced momentum equations, (8) and the expression 

for the current, (6) can be solved for the electric field. It is 
standard practice to write the result in the form 

(J -6)6 JAB (JAB)AB,,^, 
E = -qAB + ro^-g^+ri-g--r2^ ^ (14) 



3.4 Numerical Stability 

The presence of the diffusive terms on the right hand side 
of equation (15) has some important implications for any 
numerical scheme, the most obvious of which is that the 
stable time step for any explicit scheme will be proportional 
to the square of the mesh spacing. What is less obvious is 
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that the HaU term places an even more severe restriction on 
the stable time step for explicit schemes. 

In order to see this, suppose that the diffusive flux in 
equation (15) is much larger than the hyperbolic flux. Since 
ro is much smaller than the other resistivities provided that 
there is a significant charge density of particles with large (3, 
we can also ignore the term roRo in the resistance matrix. 

If we now linearise (15) about a state in which = 
and the field makes an angle 9 with the x-axis, then 



R = 



r2 



— ri cos f 



ri cos 6 



r2 cos 



Now let be a numerical approximation to 

B(t„,j'Aa;) and consider the obvious scheme 

1 

At"-"' 



^(37+1 -B7) 



^R[(l-/x)(B?+i-2B,"+B7_i) 



(16) 



where < /i < 1. This scheme is purely explicit for jj, — 0, 
purely implicit for = 1 and time centred (Crank-Nicolson) 
for fj, = 1/2. For a numerical wave of the form 

B" = B"-exp{iujj), 

we get 

B"+i = AB", 

where the amplification matrix, A is 

A = [I + ui^-R]~^[I - v{l - m)R], 

with 
_ 2At(l - cosw) 

and I the 2x2 identity matrix. Clearly the most restrictive 
condition on the timestep is for the ±1 mode, cosa; = — 1, 
for which 



4At 



(17) 



3.4.-1 Explicit Scheme (fj, = 0) 
In this case the eigenvalues of A are given by 
+ {i^r2 cos^e- 2 + 1/7-2) A 

+ 1 — !/r2 cos^ 9 — 1^2 + v^r\^ cos^ 9 + i/^r2^ cos^ 9 = Q. 

This has real roots if 

1 2 
ri cos 9 < -r2 sin 9, 

in which case stability requires 
4 



1/ < 



[r2 cos2 9 + r2 + V{r1 sin* 9 - 4rf cos^ 61)] ' 

This increases with ri cos 0, so that the Hall term in- 
creases the stable time step, but when the eigenvalues of A 
become complex, the stability condition is 

r2(l +cos'(9) 



V < 



This tells us that the stable time step tends to zero as the 
Hall term becomes large compared to ambi-polar diffusion, 
which means that an explicit scheme is very inefficient in 
such cases. This would seem to be the explanation for the 
severe restrictions on the stable timestep oxporicncod by 
HoUerbach & Riidiger (2002) in their calculations of the Hall 
effect in neutron stars. 



3.4-2 Crank-Nicolson Scheme (fj, = 1/2) 

The eigenvalues of A arc now given by 

(2i/r2 cos^ 9 + 2ur2 + 4 + i/^rf cos^ 9 + y'^rl cos^ 9)\^ 
+ (2i/Vf cos^ 9 + 2v^rl cos^ 9 - 8)A 
+ 4 + z/^rf cos^ 9 - 2vr2 cos^ 9 - 2^2 + v'^rl cos^ 9 ■ 



whose roots are 
A = 



+ ri) ± 2v^{rl sin"* 9 - 4ri cos^ 9) 



4 + 2i/r2(l 



+ cos^ 9{r'f + rj) 



As we would expect, |A| < 1 for all v, so the scheme is 
unconditionally stable. 

One might wonder whether there is some other explicit 
scheme with better stability properties. That this is unlikely 
can be seen by looking at the nature of the equation when 
the Hall term is dominant. We then get a dispersive wave 
equation for whistler waves with 

^'='^'^\ 

so that both the phase and group velocities tend to infinity 
as the wavenumber tends to infinity. This not only explains 
why the first order explicit approximation, (16), to the Hall 
term is unstable, but also suggests that there is no stable ex- 
plicit approximation. For example, the second order explicit 
scheme described by Sano & Stone (2002) is unconditionally 
unstable when the Hall term is dominant, although this was 
not apparent in their test calculation because the numerical 
resolution was low enough for the scheme to be stabilized 
by the hyperbolic term. 

Ciolek & Roberge (2002) do not use (15) to advance the 
field, instead they write E in terms of the velocity, qs of the 



E = -qa A B, 

which is valid since jS^ 1. They then use the Faraday 
equation, (4), to obtain an evolution equation for B. This 
has the disadvantage that, since there is no simple way of 
making such a scheme implicit, it requires a very small time 
step at high resolution, particularly when the Hall term is 
large compared to ambi-polar diffusion. 



3.5 Numerical Scheme for the Charged Fluids 
and Fields 

In order to obtain the state of the charged fluids and flelds, 
V„+i/2j at the half time, we flrst advance the magnetic fleld 
with the flrst order scheme 



1 .„n+l/2 _ ^ ^(M7+i/2 - MU/2) = 



At 



(rf + r|) cos^ 9 ' 



(18) 



2B 



n+1/2 



+ B 



n+1/2. 



© 0000 RAS, MNRAS 000, 000-000 



6 S. A. E. G. Falle 



where R" = R(V"). Note that the term on the right hand 
side are calculated implicitly, whereas the hyperbolic term 
is calculated explicitly. This is a block tridiagonal equation 
for B""'"'^/^ with the blocks consisting of 2 x 2 matrices and 
can readily be solved by Gaussian elimination. 

It would be nice to use a Riemann problem to compute 
the hyperbolic flux, M, but this is not possible because we 
would have to solve a Riemann problem in which the mag- 
netic field does not exert a force on the gas. In general the 
solution to such problems contains discontirmitics in the tan- 
gential velocities which are incompatible with the induction 
equation. We therefore have to be content with a centred 
approximation to the hyperbolic flux 

M7+i/2 = ^(M7+i+M,"). 

This is perfectly satisfactory as long as the resistive terms 
and numerical resolution arc such as that the magnetic field 
appears continuous on the grid. 

The densities of the charged fluids can be calculated 
from an explicit upwind approximation to the continuity 
equations and the current from a centred approximation to 
the Ampere's law, (5), using the field at the half time. Given 
J, B and the state of the neutral fluid at the half time, we 
can calculate y"^^^'^ from the reduced momentum equa- 
tions, (8), the expression for the current, (6), and the re- 
duced energy equations, (9). Note that these equations have 
to be solved by iteration if the interaction coefficients, Ku, 
depend upon the velocities, but since the solution at the old 
time provides a very good initial guess, this is not expensive. 

V"^^^^ can then be used to calculate the source term 
at the half time so that the neutral solution can be advanced 
to the full time using (13) with the second order fluxes. The 
densities of the charged fluids can also be advanced using 
■yn+\/2 ^ second order approximation to their continuity 
equations. As we shall see, this algorithm is well behaved 
even when some of the charged species have very large Hall 
parameters. 

The magnetic field is advanced explicitly to the full time 
using fluxes computed from the solution at the half time 

At ^ ^ ' Ax ^ J + l/2 3-1/2 I 

1 pn-H/2^pn-H/2 „TDn-|-l/2 „n-|-l/2^ 

Although this is explicit, the fact that the field at the half 
time has been calculated implicitly ensures the same stabil- 
ity properties as a Crank-Nicolson. Finally the neutral solu- 
tion, charged fluid densities and magnetic fleld can be used 
to calculate the charged fluid velocities and temperatures at 
the full time in the same way as for V""*"^^^. 



4 TEST CALCULATIONS 

In order to test the code described in the previous section, 
we compare the results for steady shocks with those obtained 
by a direct solution of the steady equations. The numerical 
solution was obtained by imposing the upstream and down- 
stream states at the right and left boundaries of the domain 
with a discontinuity in the middle of the domain and then 
integrating until the steady state is reached. 



4.1 Steady Solutions 

In order to calculate the steady solution, we start by trans- 
forming the induction equation, (15), to a frame in which 
the shock is steady and then setting the time derivative to 
zero. This gives 

dM _ d_^dS 

dx dx dx ' 

which can bo integrated to give 

M-Ml =M-Mfl = R— , (19) 

dx 

where the suffices L, R, denote the upstream and down- 
stream states respectively. The shock relations ensure that 
M_t ~ Mjj so that the upstream and downstream states are 
fixed points of this equation. 

We confine ourselves to the case where the neutral gas is 
isothermal, so that the steady versions of the neutral equa- 
tions can be all be integrated to give 

pu = Q, pw^ + a^p -\- iS^ = 

PUV — B^By = Py, PUW — B^B^ = Pj 

where Q, Px, Py and Pz are constants and a is the isother- 
mal neutral sound speed. Here we have used Ampere's law, 
(5), to write J in terms of B in the momentum equation. 
These equations can be solved to give p and q = {u, v, w) as 
functions of B. 

The charged fluid densities are given by 

PiUi = Qi, (20) 

where the Qi are constants. The reduced momentum equa- 
tions for the charged fluids, (8), give us 3N — 3 equations for 

the 3N — 3 components of the charged fluid velocities and 
the 3 components of the electric field. In the steady case, the 
Faraday equation, (4), tells us that Ey and are constant 
and equal to their values in the upstream state, so that only 
Ex is unknown. We therefore have 3N — 2 unknowns and 
3A'^ — 3 equations. The remaining equation can be obtained 
from (20) and the charge neutrality condition, (7). 

Given B, we can therefore calculate the neutral veloc- 
ities and the charged densities and hence M and the re- 
sistance matrix, R. We shall confine ourselves to cases for 
which the dowustrcani state is a saddle point and the up- 
stream state is a sink, in which case the steady solution 
can be obtained by integrating (19) from downstream to 
upstream. Since the gas is isothermal, any sub-shock must 
be at the downstream end of the structure. It is therefore 
a simple matter to insert a subshock and then to integrate 
from its upstream side towards the upstream state (see case 
C). 

4.2 Case A (Negligible Hall EflFect) 

This has two charged fiuids, both of which have /3 ^ 1, 
so that the Hall resistivity is neglible. The parameters are 
given in table 1. This corresponds to an oblique fast shock 
with an upstream Mach number relative to the fast speed 
Mf = 1.5 and a very small neutral pressure. The charge 
to mass ratios, Qi, and collision coefficients, Ki^i are such 
that P2 = -5.831 10^ P3 = -5.831 10^ in the upstream 
state, which is appropriate for electrons and ions in material 
with UH = 10® and B = 10"^ G (Wardle 1998). This gives 
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Table 1. Parameters for the Test Calculations 

Case A 



Upstream 


P - 


= 1 




q 


(-1.7510, 


0,0) 


B = 


(1,0.6,0) 


P2 


5 10"* 


Pz 


= 10-3 


Downstream 


P = 


= 1. 


7942 


'1 = 


(-0.9759, 


-0.6561,0) 


B = 


(1, 1.74885,0) 


, = 


8.9712 10"* 


P3 


= 1.7942 10 




a2 




-2 IOI2 


03 = 


= 108 


K12 


= 4 10^ 




= 2 10* 


a = 


0.1 


Case B 


























Upstream 


as 


for 


case A 




















Downstream 


as 


for 


case A 






















0:2 




-2 10^ 


as = 


= 10^ 




K12 


= 4 10^ 


Kl3 


= 5 105 


a = 


0.1 


Case C 


























Upstream 


P = 


= 1 




q = 


(-6.7202, 


0,0) 


B = 


(1,0.6,0) 


P2 = 


5 10-8 


P3 


= 10-3 


Downstream 


P = 


= 10.421 


q = 


(-0.6449, 


-1.0934, 0) 


B = 


(1,7.9481,0) 


P2 = 


5.2104 10"'^ 


pa 


= 1.0421 10 




a2 




-2 I012 


03 = 


= 108 


i^l2 


= 4 10^ 




= 2 10-* 


a = 


1 




0.0 



X 



0.2 



0.4 



0.6 



0.8 



X 



Figure 1. Neutral x velocity and y component of magnetic field 
for case A with Ax = 5 10^3. The line is the solution to the 
steady equations and the markers are the time dependent numer- 
ical solution. 



Ui 
-1.0 

-1.2 

-1.4 

-1.6 



0.0 



0.2 



0.4 



0.6 



0.8 



X 



Figure 2. Neutral x velocity for case A with Ax = 2.5 IQ-^. The 
line is the solution to the steady equations and the markers are 
the time dependent numerical solution. 



ro = 2 10"^^, ri = 1.16 10"''' and ra = 0.068, so that the 
Hall term is negligible compared with ambi-polar diffusion. 

EYom Fig 1, which shows the neutral x velocity and 
the y component of the magnetic field, it is clear that the 
agreement between the two solutions is excellent at high 
resolution. Even at the much lower resolution shown in Fig 
2, the errors are very small. 

In order to make this quantitative, we define the Li 
error in a primitive variable, P, by 



1 

in - Jo) ^ ' 



(21) 



where Pj is the numerical solution in cell j and Ps{x) is 
the solution calculated from the steady equations. Here jo 
and ji are such that the region xq ~ jo Ax to xi — jiAx 
covers the shock structure, but does not include much of 
the upstream and downstream uniform regions. Since the 
numerical shock position depends upon the evolution from 
the initial data, we impose a shift in x on the steady solution 
so as to minimize the error. 

For case A, with xo = —0.06, xi = 0.94, the error in 
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X 

Figure 3. Neutral x velocity and y,z components of magnetic 
field for case B Ax = 2 10~^. The line is the solution to the 
steady equations and the markers are the time dependent numer- 
ical solution. 




X 

Figure 4. Neutral x velocity for case B with Aa; = 5 IQ-^. The 
line is the solution to the steady equations and the markers are 
the time dependent numerical solution. 



the neutral x velocity is e„j = 5.68 10^ ' for Ax = 0.005 
and e„i = 1.26 10~^ for Ax = 0.025. This corresponds to 
eui oc Aa;^'^^, i.e. very nearly second order convergence. For 
Ax = 0.0125, the corresponding error is 3.52 10"*, which, 
when compared to Ax — 0.005, gives a convergence rate 
of tui oc Ax^'^^. This indicates that the asymptotic con- 
vergence is indeed second order, as we would expect for a 
smooth solution. 



4.3 Case B (Large Hall Effect) 

This case also has two charged fluids with p2 = —5.831 10® 
as in case A, but with Ps = 0.2332. This gives ro = 2 10"^, 
ri = 0.0116 and r2 = 0.00272, so that the Hall term dom- 
inates the ambi-polar term. In all other respects, the up- 
stream and downstream states are for Case A. 

This calculation is more demanding than case A since 
the fact that the Hall term is dispersive means that the 
upstream state is now a stable spiral instead of a stable 
node. The shock structure therefore contains large oscilla- 
tions whose wavelength is significantly smaller than the total 
width of the shock structure. Nevertheless, Fig 3 shows that 
at sufHciently high resolution, the agreement between the 
solutions is excellent and Fig 4 shows that it is still satis- 
factory even when there axe only ~ 50 mesh points in the 
shock structure. Note that, as we would expect, there is a 
significant z component of the magnetic field in this case. 

In this case, with xo = 0, xi = 0.4, the error in the 
neutral x velocity is = 3.79 lO^** for Ax = 0.002 and 

= 2.37 10"^ for Ax = 0.005. This gives e^^ oc Ax^ °, 
which shows that at we have already reached second order 
convergence rate at these resolutions. Again, this is what 
we would expect for a smooth solution. Note, however, that 
since the solution contains more structure than case A, the 
errors are larger even though the resolution is higher. 



© 0000 RAS, MNRAS 000, 000-000 



Multi-fluid Magnetohydrodynamics 9 




0.05 



0.05 



0.10 



0.15 
X 



0.20 




0.10 



0.15 
X 



0.20 




0.05 



Figure 5. Neutral x velocity, y component of magnetic field and 
fluid 2 X velocity for case C with Aa; = 10~^. The line is the 

solution to the steady equations and the markers are the time 
dependent numerical solution. 




Figure 6. Neutral x velocity for case C with Aa; = 5 10~^. The 
line is the solution to the steady equations and the markers are 
the time dependent numerical solution. 



4.4 Case C (Neutral Sub-shock) 

This is similar to case A except that the neutral sound speed 
a = 1 and the upstream fast Mach number M/ = 5, which 
mean that there is a neutral sub-shock at the downstream 

end of the structure. Fig 5 shows that the agreement between 
the two solutions is again excellent, except for the finite 
width of the sub-shock in the solution computed with the 
time dependent code. There is also an error in the charged 
fluid X velocities inside the sub-shock, but this does not af- 
fect the solution elsewhere. This error arises because of the 
exact steady solution has a discontinuity in the electric field 
at the sub-shock and could be reduced by adding viscosity 
to increase the width of the sub-shock. However, there is ht- 
tle point in doing this since the assumption that the inertia 
of the charged particles is negligible is not valid within the 
sub-shock. Note that in this case also, the results are still 
satisfactory at much lower resolution (Fig 6). 

Since the neutral x- velocity is discontinuous at the neu- 
tral sub-shock, which is smeared out in the numerical solu- 
tion, we expect the rate of convergence to be first order. In 
fact, with xo = 0.065, xi = 0.215, we get e„i = 6.88 10"^ for 
Aa; = 0.001 and = 6.98 10"^ for Ax = 0.005. This gives 
oc Ax''"^^, which is sigruficantly better than first order. 
This merely indicates that, at this resolution, the error is not 
entirely dominated by that at the sub-shock. If we increase 
the resolution to Ax = 5 10~*, we get = 3.94 10~^ which 
corresponds to Sui oc Ax°'^^ when compared to Ax = 0.001. 
That this is somewhat worse than first order is probably be- 
cause rounding errors are becoming significant at such high 
resolution. 



5 CONCLUSIONS 

It is clear that the code described in this paper is both ro- 
bust and accurate and it also has the advantage that, unlike 
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explicit methods, the timestep at high numerical resolution 
is not restricted by either the Hall or arnbi-polar terms. Al- 
though this is particularly important when the Hall term 
dominates over ambi-polar diffusion, it also gives a signifi- 
cant increase in efficiency at high resolution when this is not 
the case. A code that is efficient at high resolution is neces- 
sary if one wishes to study multi-fluid shock structures with 
a distribution of grain sizes and realistic physics. The fact 
that a Hall term that is large compared to ambi-polar dif- 
fusion imposes such a severe restriction on the timestep for 
explicit schemes also has important implications for any nu- 
merical calculations that include this effect. It may be that 
at relatively low resolution an explicit scheme is stabilized 
by the rmmcrical dissipation due to the hyperbolic terms, 
but such a scheme will require a very small timestep if the 
resolution is sufficient to give accurate results. 

Although we have only described a one dimensional al- 
gorithm, its structure is such that is a simple matter to 
extend it to multi-dimensions. A multi-dimensional version 
of this code would have numerous applications to the dy- 
namics of molecular clouds in general and star formation in 
particular. For example, Wardle (1991) has suggested that 
the instability that affects C-typc shocks is generic and is 
likely be important whenever there there is both ambi-polar 
diffusion and a dynamically significant magnetic field. It can 
also be applied to accretion discs, such as proto-planetary 
discs and those in dwarf novae, in which the high density re- 
duces the ion Hall parameter to the point where both ambi- 
polar diffusion and the Hall effect are significant (e.g. Sano 
& Stone 2002). 

We have not considered the other source terms such 
as the ionization and recombination, chemical reactions and 
radiative cooling, but these terms do not present any obvious 
difficulty. There may be occasions when they are stiff, but 
this merely requires locally implicit methods rather than the 
globally implicit method that we have used for the induction 
equation. 
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